# Write your code here
load(file = url("https://github.com/StatisticsSU/SDA1/blob/main/assignments/assignment1/Boston_census_data.RData?raw=true")) Statistik och dataanalys I, 15 hp
Inlämningsuppgift 1
Warning
Den här inlämningsuppgiften förutsätter att följande paket finns installerade:
mosaicdplyrgeosphereleaflet
Paket kan installeras via kommandot install.packages('packagename'), där 'packagename' är namnet på paketet, t.ex 'mosaic'.
Introduktion
I den första inlämningsuppgiften ska ni självständigt i grupper om tre analysera ett dataset i programmeringsspråket R. Till skillnad från datorlaborationerna finns det minimalt med kodexempel. Datorlaborationerna går igenom de flesta momenten som behandlas i inlämningsuppgiften, så se till att göra klart dessa innan.
Instruktioner
I denna inlämningsuppgift ska ni analysera ett datamaterial som beskriver ca 500 distrikt i Boston år 1970. Datasetet är en modifierad version1 av originaldata2 som användes i en studie3 där författarna predikterade medianhuspriser i olika distrikt med hjälp av en rad förklaringsvariabler.
Följande variabler finns i datasetet boston_census_data.Rdata (ladda ner), som innehåller 480 observationer. Notera att en observation motsvarar ett distrikt:
town: Stadsdel.longitude: Longitud koordinat.latitude: Latitud koordinat.median_home_value: Medianhuspriset (enhet 1K USD).crime_rate: Brott (per 1000 invånare).zoned_25k_p: Andel av stadsdelens bostadsmark ämnad för marklotter större än 25000 kvadratfot.indust_p: Andel tunnland ägd av företag utanför detaljhandel.borders_charles: Charles River dummy variabel (= 1 om området gränsar till floden, 0 annars).NOx: Koncentration av kväveoxider (andelar per 10 miljon).n_rooms_avg: Genomsnitt antal rum i ägda bostäder.before_1940_p: Andel ägda bostäder byggda före 1940.employ_dist: Viktat avstånd till fem arbetsförmedlingscentra i Boston.radial_access: Index som mäter tillgång till stadsmotorvägar.tax_rate: Fastighetsskatt per 10000 USD.pupil_teacher_ratio: Lärartäthet mätt som elev per lärare.lower_stat_pct: Procentandel med låg socioekonomisk status i termer av utbildning eller arbete.dist_fenway_park: Avstånd till stadion Fenway Park.
Inlämningsuppgiften ska lämnas in i form av ett html-dokument genererat av Quarto. Kontrollera att ni inte får några felmeddelande när du genererar HTML-dokumentet. Läs sedan igenom HTML-dokumentet noggrant innan ni lämnar in det. Använd tydliga figurer och namnge axlarna med tydliga variabelnamn. Glöm inte att skriva era namn i början av dokumentet, där det nu står Namn 1, Namn 2 och Namn 3.
0. Ladda in data
💪 Uppgift 0.1
Ladda in datasetet Boston_census_data med följande kod.
1. Kriminalitet i Boston
I detta avsnitt ska ni analysera kriminaliteten i Boston med hjälp av variabeln crime_rate.
💪 Uppgift 1.1
Vad kan man generellt säga om kriminaliteten i Boston 1970? Använd lämpliga figurer och mått för att ge en beskrivning.
Uppgift 1.1 - Svar
Variabeln crime_rate mäter brott per 1000 invånare. Vi börjar med att använda kommandot “favstats” för att få en bra överblick över variabeln. Favstats säger oss att det minsta värdet är 0.00632 och 88.9762 är det högsta, det ger en varians på 88,96988 vilket i sammanhanget är väldigt mycket och säger oss att brottsligheten skiljer sig väldigt mycket åt mellan olika distrikt. Medelvärdet är 3.66526 med en standardavvikelse på 8.746 medan Medianen enbart är 0.253715 med en IQR på 3.599. I och med att datasetet är snedvridet så anser vi att medianen är ett bättre mått för denna variabel.
Genom att använda en boxplot bekräftar det tanken på att datasetet är väldigt snedvridet då vi ser väldigt många outliers. Vi kan även se om vi gör ett histogram istället att majoriteten av distrikten, över 80%, har en crime_rate på mellan 0 och 10 brott per tusen invånare, vilket med andra ord betyder att de flesta distrikten har en relativt låg crime_rate medan några få distrikt har en väldigt hög crime_rate.
# Write your code here
suppressMessages(library(mosaic))
suppressMessages(library(dplyr))
suppressMessages(library(geosphere))
suppressMessages(library(leaflet))
favstats(~crime_rate, data = Boston_census_data) min Q1 median Q3 max mean sd n missing
0.00632 0.0829725 0.253715 3.681942 88.9762 3.66526 8.746156 480 0bwplot(~crime_rate, data=Boston_census_data)
histogram(~crime_rate, data = Boston_census_data, breaks = 10)
💪 Uppgift 1.2
Distrikten tillhör olika stadsdelar som anges i den kategoriska variabeln town? Det finns 88 olika sådana stadsdelar (towns).
Skiljer sig brottsligheten åt mellan de olika stadsdelarna? Undersök stadsdelarna Boston East Boston, Boston Downtown, Cambridge, samt ytterligare en stadsdel som ni själva väljer. Besvara frågan med hjälp av lämpligt valda figurer och statistiska mått.
Tip
Innan ni påbörjar analysen, skapa en ny data frame som enbart innehåller de stadsdelar som ni vill jämföra. Det kan göras till exempel med funktionen filter().
Uppgift 1.2 - Svar
Efter att vi har skapat vårt filter som enbart innehåller de utvalda stadsdelarna, n=63 istället för 480, kan precis som på förra deluppgiften börja med kommandot favstats. Där ser vi att variansen fortfarande är väldigt stor 51.1358 - 0.29819 = 50.83761 samt att medianen och medelvärdet skiljer sig åt, median = 2.3139 och medelvärdet = 6.261704. Det säger oss att fördelningen fortfarande är skev, både medelvärdet och medianen är högre här än om vi tittar på alla observationer.
Tittar vi på histogrammet ser vi tydligt att den är högerskev med en markant outlier med en crime_rate på 50 per 1000 invånare. Fördelningen är unimodal med ett typvärde på en crime_rate mellan 0 och 10 per 1000 invånare.
I boxploten ser vi att crime_raten är minst i Newton med väldigt liten spridning medan den störst i Boston Downtown med en median på runt 17-18, det finns också en outliner på över 50.
# Write your code here
Boston_census_new <- filter(Boston_census_data, town == "Boston East Boston" |town == "Boston Downtown" |town == "Cambridge" |town == "Newton")
favstats(~crime_rate, data = Boston_census_new) min Q1 median Q3 max mean sd n missing
0.29819 0.61913 2.3139 8.07211 51.1358 6.261704 9.129887 63 0histogram(~crime_rate, data = Boston_census_new)
bwplot(crime_rate ~ town, data = Boston_census_new)
💪 Uppgift 1.3
Vilka två variabler i datasetet Boston_census_data korrelerar mest med brottslighet? Beskriv sambandet mellan brottslighet och var och en av dessa två variabler.
Tip
Kom ihåg att korrelation mäter det linjära sambandet mellan numeriska variabler.
Uppgift 1.3 - Svar
Vi skapar ett nytt dataset som enbart innehåller numeriska variabler som vi döper till “Boston_census_data_15_variables”. Sen för att hitta vilka två variabler som korrelerar mest med brottslighet gör vi en correlation matrix, denna är dock svåröversiktlig så då är det bättre att göra en corrplot, där man tydligare ser vad som korrelerar med vad. Detta kan man senare dubbelkolla med correlation matrixen, ifall vissa pluppar är snarlika varandra.
Vi hittar 3 variabler som ser ut att korrelera mest med crime_rate: radial_access, median_home_value och tax_rate. Vi använder funktionen cor() för att se vilka två som har högst samband. radial_access: 0.62, median_home_value: -0.45 och tax_rate: 0.58.
Radial_access och crime_rate har ett relativt starkt positivt samband på 0.62, det innebär ju större tillgång till stadsmotorvägar desto högre crime_rate.
tax_rate och crime_rate har också ett postivt samband på 0.58, det innebär att ju högre tax_rate desto mer brott begås per 1000 invånare.
# Write your code here
Boston_census_data_15_variables <- Boston_census_data[, c("longitude","latitude","median_home_value","crime_rate","zoned_25k_p","indust_p","NOx","n_rooms_avg","before_1940_p","employ_dist","radial_access","tax_rate","pupil_teacher_ratio","lower_stat_pct","dist_fenway_park")]
correlation_matrix_Boston <- cor(Boston_census_data_15_variables)
round(correlation_matrix_Boston, 3) longitude latitude median_home_value crime_rate zoned_25k_p
longitude 1.000 0.140 -0.326 0.056 -0.168
latitude 0.140 1.000 0.015 -0.092 -0.123
median_home_value -0.326 0.015 1.000 -0.450 0.395
crime_rate 0.056 -0.092 -0.450 1.000 -0.196
zoned_25k_p -0.168 -0.123 0.395 -0.196 1.000
indust_p 0.035 -0.063 -0.600 0.403 -0.522
NOx 0.138 -0.085 -0.524 0.415 -0.506
n_rooms_avg -0.227 -0.069 0.680 -0.213 0.299
before_1940_p 0.183 0.060 -0.488 0.350 -0.559
employ_dist 0.011 -0.069 0.368 -0.378 0.672
radial_access 0.015 -0.226 -0.479 0.623 -0.303
tax_rate 0.031 -0.189 -0.579 0.580 -0.300
pupil_teacher_ratio 0.295 -0.002 -0.515 0.283 -0.375
lower_stat_pct 0.172 0.034 -0.765 0.464 -0.424
dist_fenway_park 0.435 -0.200 -0.003 -0.123 0.400
indust_p NOx n_rooms_avg before_1940_p employ_dist
longitude 0.035 0.138 -0.227 0.183 0.011
latitude -0.063 -0.085 -0.069 0.060 -0.069
median_home_value -0.600 -0.524 0.680 -0.488 0.368
crime_rate 0.403 0.415 -0.213 0.350 -0.378
zoned_25k_p -0.522 -0.506 0.299 -0.559 0.672
indust_p 1.000 0.761 -0.407 0.632 -0.706
NOx 0.761 1.000 -0.315 0.725 -0.766
n_rooms_avg -0.407 -0.315 1.000 -0.256 0.242
before_1940_p 0.632 0.725 -0.256 1.000 -0.744
employ_dist -0.706 -0.766 0.242 -0.744 1.000
radial_access 0.592 0.608 -0.190 0.447 -0.486
tax_rate 0.717 0.666 -0.278 0.499 -0.529
pupil_teacher_ratio 0.382 0.178 -0.289 0.264 -0.237
lower_stat_pct 0.640 0.615 -0.607 0.633 -0.544
dist_fenway_park -0.394 -0.325 0.020 -0.304 0.686
radial_access tax_rate pupil_teacher_ratio lower_stat_pct
longitude 0.015 0.031 0.295 0.172
latitude -0.226 -0.189 -0.002 0.034
median_home_value -0.479 -0.579 -0.515 -0.765
crime_rate 0.623 0.580 0.283 0.464
zoned_25k_p -0.303 -0.300 -0.375 -0.424
indust_p 0.592 0.717 0.382 0.640
NOx 0.608 0.666 0.178 0.615
n_rooms_avg -0.190 -0.278 -0.289 -0.607
before_1940_p 0.447 0.499 0.264 0.633
employ_dist -0.486 -0.529 -0.237 -0.544
radial_access 1.000 0.910 0.453 0.513
tax_rate 0.910 1.000 0.455 0.572
pupil_teacher_ratio 0.453 0.455 1.000 0.364
lower_stat_pct 0.513 0.572 0.364 1.000
dist_fenway_park -0.206 -0.233 0.043 -0.161
dist_fenway_park
longitude 0.435
latitude -0.200
median_home_value -0.003
crime_rate -0.123
zoned_25k_p 0.400
indust_p -0.394
NOx -0.325
n_rooms_avg 0.020
before_1940_p -0.304
employ_dist 0.686
radial_access -0.206
tax_rate -0.233
pupil_teacher_ratio 0.043
lower_stat_pct -0.161
dist_fenway_park 1.000suppressMessages(library(corrplot))
corrplot(correlation_matrix_Boston)
cor(radial_access ~ crime_rate, data = Boston_census_data)[1] 0.623462cor(median_home_value ~ crime_rate, data = Boston_census_data)[1] -0.4496238cor(tax_rate ~ crime_rate, data = Boston_census_data)[1] 0.58024442. Fastighetsskatt i Boston
I detta avsnitt ska ni analysera fastighetsskatten i Boston med hjälp av variabeln tax_rate.
💪 Uppgift 2.1
Vad kan man generellt säga om fastighetsskatten i distrikten? Använd lämpliga figurer och mått för att beskriva fördelningen.
Uppgift 2.1 - Svar
Fastighetsskatten varierar mellan 187 och 711 per 1000 USD, vilket ger en varians på 524. Medianen och Medelvärdet skiljer sig inte lika mycket som det gjorde på brottsligheten, utan ligger här på 330 respektive 409.33. Standardavvikelsen är 168.5777 och IQR är 385.
Om vi gör ett histogram ser vi att det är en bimodal fördelning samt att det ser ut att vara två olika grupper i och med att det är ett glapp mellan staplarna. Typvärderna ligger runt 300 och 650. Antagligen finns det ett samband med en annan variabel som avgör om man betalar lägre eller högre fastighetsskatt.
# Write your code here
favstats(~tax_rate, data=Boston_census_data) min Q1 median Q3 max mean sd n missing
187 280.75 330 666 711 409.3271 168.5777 480 0histogram(~ tax_rate, data = Boston_census_data, col = "navyblue", type="count")
💪 Uppgift 2.2
Låt oss skapa en ny variabel cat_tax som anger om ett distrikt betalar låg (low), medel (medium), eller hög (high) fastighetsskatt. Vi definerar skattekategorierna enligt
low:tax_rate\(\leq\) 250,medium: 250 \(<\)tax_rate\(\leq\) 400,high:tax_rate\(>\) 400.
Följande kod skapar och lägger till variabeln cat_tax i Boston_census_data
Boston_census_data$cat_tax <- cut(Boston_census_data$tax_rate,
breaks=c(0, 250, 400, 800),
labels=c('Low', 'Medium', 'High'))Finns det ett samband mellan vilken skattekategori ett distrikt tillhör och dess angränsning till Charles River? Förklara med hjälp av lämplig tabell och figur.
Uppgift 2.2 - Svar
Som man ser i stapeldiagrammet skiljer sig skattenivåerna inte så mycket åt mellan om man gränsar till Charles River eller inte. Det man kan se är att det är mer troligt att man har men medium-tax rate om man gränsar till Charles River än om man inte gör det.
Fördelningen mellan skattenivåer om man inte gränsar till Charles River eller om man gör det kan vi se i korstabellen. Det är större sannolikhet att man inte gränsar till Charles River om man har låg (13,3% mot 10,3%) eller hög (39,69% mot 34,48%) skattekategori.
# Write your code here
tally(cat_tax ~ borders_charles, data = Boston_census_data, margins=TRUE, format="percent") borders_charles
cat_tax 0 1
Low 13.30377 10.34483
Medium 47.00665 55.17241
High 39.68958 34.48276
Total 100.00000 100.00000bargraph(~cat_tax, groups=borders_charles, data=Boston_census_data, type="percent", main = "Skattekategorier baserat på angräsning till Charles River")
💪 Uppgift 2.3
Hur många procent av alla distrikt i vår data ligger i angränsning till Charles River och tillhör en hög skattekategori? Hur stor andel av distrikten med hög skatt ligger inte i angränsning till Charles River?
💪 Uppgift 2.4
Vilka två variabler i datasetet Boston_census_data korrelerar starkast med tax_rate? Beskriv det parvisa sambandet mellan tax_rate och var och en av dessa två variabler. Vad kan vi säga om kausalitet för vart och ett av sambanden?
Tip
Kom ihåg att korrelation är ett mått på linjära samband mellan numeriska variabler.
Uppgift 2.4 - Svar
Skriv svaret här.
Vi skapar ett nytt dataset som enbart innehåller numeriska variabler som vi döper till “Boston_census_data_15_variables”. Sen för att hitta vilka två variabler som korrelerar mest med fastighetsskatt gör vi en correlation matrix, denna är dock svåröversiktlig så då är det bättre att göra en corrplot, där man tydligare ser vad som korrelerar med vad. Detta kan man senare dubbelkolla med correlation matrixen, ifall vissa pluppar är snarlika varandra.
Man ser ganska tydligt att det är två variabler som verkar korrelera mest med tax_rate och det är radial_access och indust_p. Vi använder cor()-funktionen för att se korrelationerna.
Radial_access och tax_rate har ett väldigt starkt positivt samband på 0.91. Det verkar finnas ett tydligt samband mellan större tillgång till stadsmotorvägar och distriktens fastighetsskatt.
Även indust_p och tax_rate har ett starkt positivt samband på 0.72. indust_p står för andel tunnland ägd av företag utanför detaljhandeln. Det kan man tolka att ju mer tunnland som ägs av företag utanför detaljhandeln, desto högre tax_rate.
Vi kan inte säga att någon av dessa är kausala samband, vi vet inte med säkerhet att det är t ex tillgången till stadsmotorvägar som påverkar fastighetsskatten, utan det kan alltid finnas en dold variabel utanför vårt dataset som påverkar både, och bara får det att se ut som att det finns ett samband.
# Write your code here
Boston_census_data_15_variables <- Boston_census_data[, c("longitude","latitude","median_home_value","crime_rate","zoned_25k_p","indust_p","NOx","n_rooms_avg","before_1940_p","employ_dist","radial_access","tax_rate","pupil_teacher_ratio","lower_stat_pct","dist_fenway_park")]
suppressMessages(library(corrplot))
corrplot(correlation_matrix_Boston)
cor(radial_access ~ tax_rate, data = Boston_census_data)[1] 0.9099416cor(indust_p ~ tax_rate, data = Boston_census_data)[1] 0.71710413. Avstånd till Fenway park
I detta avsnitt ska ni undersöka variabeln dist_fenway_park, som mäter avståndet mellan ett distrikt och Fenway park (stadion där basebollslaget Boston Red Sox spelar sina hemmamatcher).
Vi kan visualisera Fenway park och distrikten på en karta med hjälp av R-paketet leaflet. Följande kod visar platsen för Fenway park och distrikten för observationerna 30 och 45.
library(leaflet) # Install if not available
fenway_park_lat_long <- c(42.346462, -71.097250) # latitude and longitude for Fenway_park
Boston_map <- leaflet() %>%
addTiles() %>%
addMarkers(lat = fenway_park_lat_long[1], lng = fenway_park_lat_long[2], popup="Fenway park") %>%
addMarkers(lat = Boston_census_data$latitude[30], lng = Boston_census_data$longitude[30], popup="Observation 30") %>%
addMarkers(lat = Boston_census_data$latitude[45], lng = Boston_census_data$longitude[45], popup="Observation 45")
Boston_map # Show interactive map💪 Uppgift 3.1
Vilket distrikt i vår data har längst respektive kortast avstånd till Fenway park? Markera ut dessa distrikt i en interaktiv karta tillsammans med Fenway park.
💪 Uppgift 3.2
Finns det ett samband mellan dist_fenway_park och crime_rate?
4. Enkel linjär regression
I detta avsnitt ska ni anpassa och tolka enkla linjära regressionsmodeller.
💪 Uppgift 4.1
Anpassa en linjär regression med responsvariabeln NOx och den förklarande variabeln employ_dist. Rita den anpassade regressionslinjen tillsammans med data i en lämplig figur. Beskriv resultaten och tolka modellen. Utför en modellvalidering via en residualanalys och kommentera modellens lämplighet. Om modellen inte anses lämplig, vilka antaganden har inte varit uppfyllda?
Uppgift 4.1 - Svar
Vår responsvariabel,y, är NOx, vilket är koncentration av kväveoxider (andelar per 10 miljoner). Variabeln som vi ska använda för att försöka förklara variantionen av NOx är employ_dist, x, vilket är det viktade avståndet till 5 arbetsförmedlingscentra i Boston.
Utifrån regressionsmodellen kan vi se att interceptet är 0.7176, vilket är det värde som NOx har när employ_dist är 0. Lutningen, slopen är -0.0426 vilket behöver att varje gång employ_dist ökar med en enhet så minskar NOx med -0.0426. Detta kan man tolka som att ju längre avståndet är till arbetsförmedlingscentra desto mindre är koncentrationen av kväveoxider.
R-squared är 0.5861, vilket betyder att employ_dist kan förklara ca 59% av variationen i NOx.
Modellvalidering
Utifrån modellvalidering kan vi se att modellen inte uppfyller något av kriterierna. I residualplotten kan vi se att residulernas varians inte är konstant, dem är koncentrerade på den högra sidan av modellen samt har några outliners uppe i det vänstra hörnet.
Den är inte heller helt normalfördelad, då den inte följer det röda strecket i början eller slutet, även om de gör det på det stora hela.
# Write your code here
regressionmodel_NOx_employdist <- lm(NOx ~ employ_dist, data = Boston_census_data)
summary_model <- summary(regressionmodel_NOx_employdist)
summary_model
Call:
lm(formula = NOx ~ employ_dist, data = Boston_census_data)
Residuals:
Min 1Q Median 3Q Max
-0.12358 -0.05352 -0.01220 0.04348 0.22868
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.717614 0.007109 100.95 <2e-16 ***
employ_dist -0.042641 0.001639 -26.02 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07526 on 478 degrees of freedom
Multiple R-squared: 0.5861, Adjusted R-squared: 0.5852
F-statistic: 676.9 on 1 and 478 DF, p-value: < 2.2e-16plot(NOx ~ employ_dist, data = Boston_census_data, col = "cornflowerblue", ylim = c(0, 1))
y_hat <- predict(regressionmodel_NOx_employdist)
lines(Boston_census_data$employ_dist, y_hat, type = "p", col = "lightcoral")
abline(regressionmodel_NOx_employdist, col = "lightcoral")
legend(x = "topleft", pch = c(1, 1, NA), lty = c(NA, NA, 1), col = c("cornflowerblue", "lightcoral", "lightcoral"), legend=c("Data", "Predicted", "Fitted line"))
# residualer
resid <- residuals(regressionmodel_NOx_employdist)
head(resid) 1 2 3 4 5 6
-0.009743748 0.002582251 0.047758653 -0.039548322 0.101821788 -0.079159083 # residualanalys
Boston_census_data$res <- resid(regressionmodel_NOx_employdist)
Boston_census_data$y_hatt <- fitted(regressionmodel_NOx_employdist)
plot(Boston_census_data$res ~ Boston_census_data$y_hatt, ylab="Resid", xlab="y-hatt", main="Residplot")
abline(h=0)
qqnorm(resid, col = "cornflowerblue") # Create normal probability plot for residuals
qqline(resid, col = "red") # Add a straight line to normal probability plot 
💪 Uppgift 4.2
Använd modellen i Uppgift 4.1 för att prediktera koncentration av kväveoxider för observation 10, där employ_dist=10.5857. Beräkna vad residualen blir för denna observation.
💪 Uppgift 4.3
Transformera variablerna i Uppgift 4.1 (avgör själv vilken eller vilka av de två som behöver transformeras). Ett förslag är att använda Tukeys cirkel för att hitta lämpliga transformationer. Anpassa en ny linjär regression med de transformerade variablerna. Utför en modellvalidering (efter transformation) via en residualanalys och kommentera modellens lämplighet jämfört med modellen i Uppgift 4.1.
Uppgift 4.3 - Svar
Vi testar oss fram för att hitta den transformation som ger den rätaste linjen och kommer fram till att vi transformerar både x och y till log(x) och log(y).
Modellvalidering
Jämfört med modellvalideringen i 4.1 ser vi att variansen nu är konstant, vi kan inte se något tydligt mönster i residualplotten. Den är även mer normalfördelad än 4.1 även om residualerna följer samma mönster så är de närmare den röda linjen.
# Write your code here
plot(NOx ~ employ_dist, data = Boston_census_data,col = "cornflowerblue") #otransformerad
# Testar olika transformationer
plot(sqrt(NOx) ~ sqrt(employ_dist), data = Boston_census_data, col = "cornflowerblue")
plot(sqrt(NOx) ~ employ_dist, data = Boston_census_data, col = "cornflowerblue")
plot(NOx ~ sqrt(employ_dist), data = Boston_census_data, col = "cornflowerblue")
plot(log(NOx) ~ log(employ_dist), data = Boston_census_data, col = "pink")
plot(NOx ~ log(employ_dist), data = Boston_census_data, col = "pink")
plot(log(NOx) ~ employ_dist, data = Boston_census_data, col = "pink")
plot(sqrt(NOx) ~ log(employ_dist), data = Boston_census_data, col = "coral")
plot(log(NOx) ~ sqrt(employ_dist), data = Boston_census_data, col = "coral")
# Vi väljer att transformera både x och y till log(x) och log(y) för det ger den rätaste linjen.
lm_logNOx_vs_logemploy_dist <- lm(log(NOx) ~ log(employ_dist), data = Boston_census_data)
summary(lm_logNOx_vs_logemploy_dist)
Call:
lm(formula = log(NOx) ~ log(employ_dist), data = Boston_census_data)
Residuals:
Min 1Q Median 3Q Max
-0.19152 -0.07828 -0.01370 0.06108 0.28496
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.218492 0.011524 -18.96 <2e-16 ***
log(employ_dist) -0.327303 0.008832 -37.06 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.103 on 478 degrees of freedom
Multiple R-squared: 0.7418, Adjusted R-squared: 0.7413
F-statistic: 1373 on 1 and 478 DF, p-value: < 2.2e-16# Residualanalys
# residualer
resid <- residuals(lm_logNOx_vs_logemploy_dist)
# residualanalys
Boston_census_data$res <- resid(lm_logNOx_vs_logemploy_dist)
Boston_census_data$y_hatt <- fitted(lm_logNOx_vs_logemploy_dist)
plot(Boston_census_data$res ~ Boston_census_data$y_hatt, ylab="Resid", xlab="y-hatt", main="Residplot")
abline(h=0)
qqnorm(resid, col = "cornflowerblue") # Create normal probability plot for residuals
qqline(resid, col = "red") # Add a straight line to normal probability plot 
plot(log(NOx) ~ log(employ_dist), data = Boston_census_data, col = "cornflowerblue")
log_y_hat <- predict(lm_logNOx_vs_logemploy_dist)
orginal_skala_y <- exp(log_y_hat) #transformera tillbaka orginaldata
plot(NOx ~ employ_dist, data = Boston_census_data, col = "cornflowerblue")
lines(Boston_census_data$employ_dist, orginal_skala_y, type ="p", col = "lightcoral")
abline(lm_logNOx_vs_logemploy_dist, col = "lightcoral")
legend(x = "topleft", pch = c(0,5, 0,5, NA), lty = c(NA, NA, 0,5), col = c("cornflowerblue", "lightcoral", "lightcoral"), legend=c("Data", "Predicted", "Fitted line"))
💪 Uppgift 4.4
Plotta den anpassade regressionen från 4.3 i icke-transformerad skala tillsammans med observationerna (också i icke-transformerad skala) i en lämplig figur.
Uppgift 4.4 - Svar
Skriv svaret här.
# Write your code here
logy_hat <- predict(lm_logNOx_vs_logemploy_dist) # log scale prediction
y_hat <- exp(logy_hat) # original scale prediction
head(logy_hat) 1 2 3 4 5 6
-0.4398989 -0.7327355 -0.3583458 -0.4774478 -0.5178857 -0.7635522 head(y_hat) 1 2 3 4 5 6
0.6441015 0.4805925 0.6988314 0.6203647 0.5957789 0.4660081 plot(NOx ~ employ_dist, data = Boston_census_data, col = "cornflowerblue") # Data on original scale
lines(Boston_census_data$employ_dist, y_hat, type = "p", col = "lightcoral")
legend(x = "topleft", pch = c(1, 1), col = c("cornflowerblue", "lightcoral"), legend=c("Data", "Predicted"))
💪 Uppgift 4.5
Använd modellen i Uppgift 4.3 för att prediktera koncentration av kväveoxider i icke-transformerad skala för observation 10, där employ_dist=10.5857. Beräkna vad residualen blir för denna observation. Kommentera resultaten jämfört med Uppgift 4.2.
Tip
Tänk på att ta hänsyn till eventuella transformationer!
5. Multipel linjär regression
I detta avsnitt ska ni studera multipel linjära regression.
💪 Uppgift 5.1
Anpassa en linjär regression med responsvariabel logaritmerad median_home_value samt förklarande variabler lower_stat_pct och dummy-variabeln borders_charles. Tolka koefficienten för borders_charles.
💪 Uppgift 5.2
Ni ska nu utforma en modell som predikterar medianhuspriset median_home_value. Ni får endast använda följande förklaringsvariabler:
before_1940_pcrime_rateradial_accessNOxdist_fenway_park
Ni får själva välja hur många av variablerna som ska ingå i modellen. Ni får göra vilka transformationer ni vill av variablerna, inklusive responsvariabeln.
Pröva er fram metodiskt när ni väljer vilka variabler ni inkluderar i modellen, och när ni bestämmer vilka eventuella transformationer ni använder.
När ni utvärderar olika modeller kan ni förslagsvis börja med att jämföra adjusted R-squared. När ni med hjälp av adjusted R-squared har identifierat två eller tre modeller som ser lovande ut kan ni utvärdera dessa modeller ytterligare i ett andra steg.
I det andra steget ska ni utvärdera hur väl modellerna predikterar data som inte använts för att anpassa modellen. Ni kan välja en av två alternativa metoder:
- Dela in ert dataset i träningsdata och testdata. Anpassa modellen med hjälp av träningsdata, och utvärdera sedan på testdata. Ni kan exempelvis använda de första 350 observationerna som träningsdata och de sista 130 observationerna som testdata.
- Använd korsvalidering. Det är en något mer krävande metod, men också något bättre. Ni kan exempelvis göra korsvalidering med 4 folds (4-fold cross validation). Dela då upp ert dataset i fyra delar (del 1: observationer 1-120, del 2: observationer 121-240, del 3: observationer 241-360, del 4: observationer 361-480).
Sortera inte observationerna i Boston_census_data slumpmässigt. Ordningen är redan slumpmässig.
Tip
Tänk på att ta hänsyn till eventuell transformation av responsvariabeln. Om ni exempelvis har valt transformationen \(\log(y)\) är modellens prediktion av responsvariabeln \(\widehat{\log(y)}\). Ni måste då transformera den till \(\hat y\) i responsvariabelns originalskala med formeln \(\hat{y}=\exp\left(\widehat{\log(y)}\right)\). Sedan kan ni räkna ut residualen \(y - \hat y\).
💪 Uppgift 5.3
Gör en residualanalys av den valda modellen i Uppgift 5.2.
Uppgift 5.3 - Svar
Residualanalysen visar att residualerna är normalfördelade, de är också relativt varierade även om man kan utläsa ett visst mönster.
# Write your code here
model_2 <- lm(log(median_home_value) ~ before_1940_p + NOx, data=datatest)
model_2
Call:
lm(formula = log(median_home_value) ~ before_1940_p + NOx, data = datatest)
Coefficients:
(Intercept) before_1940_p NOx
4.025999 -0.002243 -1.485534 resid_model_2 <- residuals(model_2)
head(resid_model_2) 1 2 3 4 5 6
0.08295812 0.14293552 -0.09737965 0.34307293 0.21512746 0.04316374 plot(datatest$median_home_value, resid_model_2, xlab= "median_home_value", ylab='Residuals', col = "cornflowerblue") 
qqnorm(resid_model_2, col = "cornflowerblue") # Create normal probability plot for residuals
qqline(resid_model_2, col = "red") # Add a straight line to normal plot
💪 Uppgift 5.4
Använd modellen i Uppgift 5.2 för att prediktera medianhuspriset för observationerna i datasetet Boston_districts_to_predict (ladda ner).
Det här datasetet har endast de förklarande variablerna, dvs inte responsvariabeln. När vi rättar era inlämningsuppgifter kommer vi att jämföra era prediktioner med de faktiska medianpriserna (som vi har tillgång till).
Skriv ut dina prediktioner så att vi enkelt kan se dem när vi rättar.
Tip
Tänk på att ta hänsyn till eventuella transformationer av förklaringsvariablerna!
Footnotes
Totalundersökningen trunkerade medianhusvärdet till 50K för de censusdistrikten som låg över. Vi har tagit bort dessa censusdistrikt. Vi har också tagit bort variabler som är irrelevanta.↩︎
Harrison Jr, D., & Rubinfeld, D. L. (1978). Hedonic housing prices and the demand for clean air. Journal of Environmental Economics and Management, 5(1), 81-102.↩︎
Pace, R. K., & Gilley, O. W. (1997). Using the spatial configuration of the data to improve estimation. The Journal of Real Estate Finance and Economics, 14(3), 333-340.↩︎